Quantum Belief Propagation 
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We present an accurate numerical algorithm, called quantum belief propagation (QBP), for sim- 
ulation of one-dimensional quantum systems at non-zero temperature. The algorithm exploits the 
fact that quantum effects are short-range in these systems at non-zero temperature, decaying on a 
length scale inversely proportional to the temperature. We compare to exact results on a spin- 1/2 
Heisenberg chain. Even a very modest calculation, requiring diagonalizing only 10-by-lO matrices, 
reproduces the peak susceptibility with a relative error of less than 10~^, while more elaborate 
calculations further reduce the error. 



The fact that interactions are short-ranged in many 
physical systems is a major simplification in finding quan- 
tum ground states. The classic example of this is the 
density matrix renormalization group(DMRG)^, which 
relies on the ability to approximate the ground state 
by a matrix product state. While it has long been be- 
lieved that such an approximation is possible whenever 
there is a spectral gap, due to conformal field theory 
calculations [2| and decay of correlations 3, 4], only very 
recently has a general proof been given that such an ap- 
proximation is possible whenever there is a local Hamil- 
tonian and a gap[^. 

At non-zero temperature, the system is in a mixed 
state instead. Here, matrix product density operatorsf6) 
play the same role that matrix product states do in study- 
ing pure states, and a non-zero temperature plays a sim- 
ilar role when it comes to representing mixed states as 
matrix product density operators as does a gap for rep- 
resenting pure states as matrix product states. Indeed, it 
has been shown that good matrix product operator rep- 
resentations exist for quantum systems at any non-zero 
temperature 0. In this paper we present quantum be- 
lief propagation (QBP), another method for finding ma- 
trix product density operators for thermal states, which 
avoids the problem of Trotter error in other methods. 

Belief propagation's,] for classical systems (CBP) is es- 
sentially a Bethe-Peierls solution of a classical statistical 
mechanics model. CBP is exact on trees, and is often a 
very good approximation on lattices with few loops. At 
the end of this paper we will discuss application of the 
QBP equations to higher dimensional systems and trees, 
but for now we will focus discussion on belief propagation 
for one dimensional lattices. In this case, CBP become 
equivalent to a transfer matrix technique: one solves the 
problem on a chain of iV-sites to get a partition func- 
tion which depends on the value of the spin on the iV-th 
site. Then, one adds a coupling of the iV-th spin to one 
additional spin, traces out the A'^-th spin, arriving at a 
partition function which depends on the value of the spin 
on the A^-|- 1-st site. One proceeds in this way, iteratively 
solving longer and longer chains. 

In a quantum system, we will proceed in a similar way. 
However, now the operator coupling the A'^-th spin to 



the A^ -f 1-st spin need not commute with the rest of the 
Hamiltonian. This means that to study properties of the 
chain of A^ + 1-spins, it is not sufficient simply to know 
the density matrix of the A^-th spin in an A^-site chain. 
However, our physical intuition tells us that at a non-zero 
temperature, quantum effects should be short range. Our 
main result in the next section realizes this intuition in 
the QBP equations, which involve two terms. One is a 
"classical" term which is local, coupling the A^-th spin 
to the A^ -|- 1-st spin. The other is a "quantum" term, 
which is non-local, coupling the A^ + 1-st spin to several 
other spins; however, this quantum term is exponentially 
decaying on a length scale set by the inverse tempera- 
ture. This will allow us to accurately describe statistical 
properties of the A^ -I- 1-st spin knowing only the reduced 
density matrix of spins N,N— 1,...,N — + for some 
Iq, so that we keep track of a reduced density matrix for 
Iq spins. We then iterate these equations by tracing out 
spin A^ — ^0 + 1 and then computing statistical properties 
of spin A^ + 2, keeping always a reduced density matrix 
on the last spins on the chain. 

Quantum Belief Propagation Equations — The 
QBP equations describe how the partition function, 
exp(— /3i7), changes when H is changed by some per- 
turbation A. We take A, H Hcrmitian throughout. 
We will apply this result to the following case: we 
have a nearest-neighbor Hamiltonian, with /i^.i+i the 
term coupling spin i to spin i + 1. We will take 
i+i, so that H is the Hamiltonian 
for the first A^ spins, and A — /iat.jv+i- We define 



H + sA. 



(1) 



Then, Hi ^ h^^^ + A = h<^^+'') . 

Although ultimately we want to compute exp(— /Jif^) 
at s = 1, we begin by computing the derivative 
ds eKp(~f3Hs). Let Aab{s) denote matrix elements 
of A in a basis of eigenstates of Hg, with energies 
Ea{s), Eh{s). Define the operator ^"''^ by its matrix ele- 
ments: = Aabis)SiEa{s) - Eb{s) - io). Wc will 
do most of the calculations in terms of A'^''^ rather than 
A for simplicity, although we will convert certain results 
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back into results in terms of A itself using the integral 

(a,exp(-/3i/,)) = I dLj(^d,cxp[-P{Hs + eA^^')]y (2) 

where all derivatives with respect to e are taken at e = 
throughout. 

One result for the derivative on the right-hand side of 
Eq. ([2]) (not the result we will use!) is (^d^exp[—P{Hs + 

eA'^-'')]) = -exp(-/3i/3)(2S(M;:i^'^,'^y One prob- 
lem with this is that the operator norm || 2iP(^f£liLL^'^.s || 
may be exponentially large. The QBP equations will im- 
prove on this, writing 9^ exp(— /3i7s) — r]exp{—(3Hs) + 
e^p{~'f3Hs)ri\ where \\r)s\\ < (/3/2) HAIL Finally, rj wih 
be a local operator as discussed below|ll|. 

To find the QBP equations, we begin by studying a 
certain correlation function. Let B be an arbitrary oper- 
ator. Then, 



aetr(exp[-/3(iJ, -FeA'^'")]B 



(3) 



/ dTtr(exp(-/3i/,)A"'^(-ir,i7,)-B) 
Jo 

exp(/3w) — 1 



-tr(exp(-/3i?,)A"'''B). 

LU 

Adding and subtracting -{f3/2)tT{{exp{-PH,), A'^^'jB) 
to Eq. da]) gives 



5etr(exp[-/3(iJ, -I- eA"'")]B) 
-ftr({exp(-/3i7,),A--^}B) 



(4) 



tr(exp(-/3i73)A"^^B) 



where we used tr(A'^'* exp(— /3i7s)-B) = 
exp(/3u;)tr(exp(-/377,)A'^^"B). 

We now focus on the correlation function 
tr(exp(— /3i/s)A"'''i3), following a procedure very similar 
to that in [§]. In [§| the result was expressed in terms of 
anti-commutators, as A, B were fermionic operators and 
hence anti-commutated if they were separated in space, 
while here we will express the result in terms of commu- 
tators, since we intend to apply it to bosonic operators 
where the Lieb-Robinson bound is expressed as a bound 
on the commutator. Using tT{exp{—/3Hs)A'^'^B) = 
tr(exp(-/3i^,)[A'^'^ B]), we have 



1 — cxp(/3lj) 



1 + e 



tiiejipi-f3Hs)A'^'''B) 



/3F(u;)tr(exp(-/37^,)[A"•^B]), 



where 



1 l-f e''" 1 _ coth(/3w/2) 

2 1 - e^" + ^ - 2 



(5) 



(6) 



Thus for any operator B, we have 9ctr(exp[— -|- 
eA'^^')]B) = -ftr({exp(-/3i/«),A"'^}B) + 
/3F(w)tr(exp(-/3i7^)[A'^'",B]), and hence 

d,exp[-p{Hs + eA^^')] (7) 
- -^{exp(-/3i/,), A-^^} + /3i^H[exp(-/3H,), A-^^]. 



We define 



lis = dtj?7^ 



(8) 



and integrate Eq. ([7]) over u to get 

ds eM-l3Hs) = Vs exp{-pHs) + eM-f^H^H- (9) 

Eq. (HD is the QBP equation for dsexp{-PHs). The 
anti-commutator in Eq. ([7]) is a "classical term" . It is the 
only term present if [A, H] = 0, in which case these equa- 
tion reproduce the classical belief propagation equations. 
The commutator in Eq. ([7]) is a "quantum" term. It is 
odd in oj and vanishes at a; = 0. We now discuss local- 
ity properties of the quantum term, assuming that H is 
local in the sense of having a Lieb-Robinson bound [lot. 
We have0 f3 J dLuF{Lu)A'^--' = 



J duj(^~l^coth{(3uj/2) + ^'^A'^-' 
- I duj^{uj -27^711/(3)-^^^^' 

/oc 
At sign(t) exp(-27mi//3)A(t, Hs), 
„ . . -oo 



(10) 



where the sum is over integer n and A{t,Hs) = 
exp{iHst)Aexp{—iHst). The integral over t is exponen- 
tially decaying for t > (3 and so, using a Lieb-Robinson 
bound, T] is local. 

Numerical Implementation and Results — In this sec- 
tion we discuss the numerical implementation of the QBP 
equations, and the results of their application to the an- 
tiferromagnetic spin-1/2 Heisenberg chain. The idea is 
to take the QBP equations, which depend on A and 
H — h'-^\ and instead set H = J2f=~N-io+2 
some constant Iq > 2. Since rj is local, this approxima- 
tion is justified for small enough /3; as /? increases, Iq 
must increase and the numerical effort is exponential in 
Iq, since one must diagonalize matrices of size 2'". 

Although translation invariance is not necessary, it 
is useful as we can then apply translations to H = 
J2f=~N-io+2 -^i^i+i t° make H equal to h^^"^^\ so for 
the Heisenberg chain with coupling constant J — 1, 
H = T,^,l<^<lo-2 Si ■ Si+i. We set A = hi^-i,i„. We 
set T to be the operator which translates one site to the 
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right. We define p to be a reduced density matrix for 
sites l...lo, so that p is a 2'" dimensional matrix. Define 

= 5'cxp(/ f],,ds'), (11) 

where the S' denotes that the integral is s'-ordered. 

The algorithm to compute the free energy per site of an 
infinite chain proceeds through the following four steps. 
(1) Initialize p to exp{—/3H). (2) Approximately com- 
pute the operator O as described below. (3) Go through 
a series of rin iterations of the following three steps: 

(A) replace p with OpO^ . (B) trace out the first site, so 
that p is replaced with tri(p) ® i-ig+i- Here, tri denotes 
a partial trace over the first site, and lip+i is the unit 
operator on the Iq + 1-th site. At this stage p is now a 
density operator on sites 2...lo + 1. Translate by one site, 
so that p becomes a density operator on sites I...I0 again. 
(C) Define Z = ti{p) and then replace p with Z~^p. 

After these iterations, we (4) Output ln(Z) from the 
last step as the free energy. This procedure relies on a 
series of iterations of steps (A), (B), (C) to find the p 
which is the fixed point of the map 

p^ Z-^T^[tYi{OpO^)<E,li,+iy. (12) 

The number of iterations required for convergence ap- 
pears to increase roughly linearly with Iq and /3, but the 
computational effort grows only linearly in the number 
of iterations. The locality of the operators rj, O justifies 
tracing out the first site in step (B). After riit iterations, p 
is approximately proportional to the reduced density ma- 
trix tri „;j(exp(— /3/i('''~^+"'*))), where the partial trace 
is over sites L.-nn. 

To compute correlation functions of operators such as 
SfS^^j we follow the following procedure. We go through 
steps (1), (2), (3) as above, with nu iterations in step (3). 
On the last of the na iterations, after step (A) we copy 
p to a new matrix Pcorr- We then replace Pcorr with 
Sgpcorr, thus inserting the first of the two operators. On 
steps (B,C) of the last iteration we replace of the map 
p^ Z~^tri{p) ^1 and pcorr ^~^tri(pcorr) <8> 1, where 
the same Z is used for both p and Pcorr- When then pro- 
ceed through j more iterations of steps (A),(B),(C), and 
on the last iteration after step (A) we replace Pcorr with 
Sq Pcorr, inserting the second of the two operators, before 
proceeding with steps (B),(C). We then proceed through 
several more iterations in which at each step we map 

p Z-^ili{OpO^) ® 1, Pcorr Z~^tTi{OpcorrO^) ® 1 

and finally output the ratio tr(pcorr)/tr(p). 

To compute the matrix O, we approximate by dividing 
the integral over s' into UsUce different slices: 

C' « exp( — ^ ^)...exp( — ^)exp( — —j, (13) 

'^slice ^slice '^slice 

where s(m) = (m — 1/2) /usUce- To compute 
the matrix exponential for each slice, exp(r/s/) for 



S' = (l/2)(l/n«Kee),(3/2)(l/n,i,ce),(5/2)(l/n«Hce),... 

we used a Taylor series method, while to compute rjs' 
itself we diagonalize Hg' and transform A into a basis of 
eigenvector of Hgi. We use the fact that H and A con- 
serve total 5*^ to speed up both this diagonalization and 
the multiplication p — > OpO^ . 

Since rjs is non-Hermitian, but WrjsW is not too large, 
the Taylor series method is a good choice for computing 
the matrix exponential. Eq. (jl3|) approximates the s'- 
dependence of rjs' by taking its value halfway through the 
slice, which gives us first order accuracy in ds'rjs> for free. 
Further, rjs' is in fact only very weakly dependent on s' 
and so a small risUce suffices, as seen by the following test: 
set Iq = 3. Then, after one iteration of steps (A),(B) and 
before normalizing in step (C), the matrix p should be 
equal to, the thermal density matrix for a three site chain. 
At (3 = 1, the largest eigenvalue should be equal to e. For 
a calculation with risUce — 1, the largest eigenvalue was 
found to be 2.718224; for risUce = 2, we find 2.718267, 
and for ngUce = 3 we find 2.718275. 

We have tested QBP by computing the susceptibil- 
ity, 0y^-j( Sj+ i sX the susceptibility peak, using the 
known location|12l| of the peak at Tmax = = 
0.64085103085. The exact resuh for the susceptibility 
is 

\exact (T^ax) = 0.146926279... (14) 

while a calculation using Iq — 5, nn — 20, UgUcc — 20, 
and correlations up to j = 20 gives 

XQBp{Trnax) = 0.146927031... (15) 

for a relative error of « 5 * 10^^. The calculation took 
< 0.08 seconds on a 1.5 GHz PowerPC G4 processor, 
and using conservation of the largest matrix diago- 
nalized was 10-by-lO. A larger calculation, with Iq = 9, 
nsUce = 50, riit = 30 improves this to XQBp{Tmax) = 
0.146926251... for a relative error of « 2 * 10"'^. 

We calculated the specific heat C as a function of tem- 
perature in two ways: first, by calculating /3^9^ In(^) us- 
ing ln(Z) from the algorithm, and, second, by calculating 
—?>(5'^dp{S^Sf^i). The results are shown in Fig. 1, where 
to take derivatives we calculated ln(Z) and {SfSf^^) for 
/3 = 0.1, 0.2, 0.3, 10.0. As Iq gets larger, the curves re- 
main accurate to lower temperature. The peak specific 
heat for Iq = 1 was 0.349914... from the second derivative 
calculation and 0.349717... from the first derivative, both 
of which compare very well with the Bethe ansatz result 
of 0.349712.... 

The accuracy can be improved by going to larger 
Iq. Another improvement is to take H — ft,('o-2) _|_ 

(l/2)/i;„_2,io-i and A = [l/2){hi^^2,ia~i + ^io-i.'o) i^i- 
stcad oi H = /i^'"^^) and A — hig-i^ig. We still have 
H + A = T^HT + hQ^i in this case, but the slightly differ- 
ent form of the perturbation seems to work better. The 
figure inset shows a comparison of Bethe ansatz data to 
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Temperature 



FIG. 1: Specific lieat against temperature for lo ~ S (dasfied 
line), 5 (dotted line), 7 (solid line). Curves that go negative 
are from —3f3^dp{SfSi^i), while those that diverge positively 
are from l3^dpln{Z). Inset: lo ~ 9 and Bethe ansatz. 



lo = 9 (where the largest matrix diagonalized is 126 di- 
mensional). 

Discussion — The implementation of the QBP equa- 
tions here must be considered as preliminary. More work 
is needed to optimize the algorithm and, most impor- 
tantly, to quantify the sources of error. Despite this, 
the method yields accurate results, giving qualitatively 
correct behavior even at Zq = 3 where QBP can be im- 
plemented as a "pen-and-paper" technique. 

In contrast to QBP, thermodynamic DMRGflsj com- 
putes low-lying eigenvalues for finite size chains up to size 
N = 30; the results were accurate to quite low temper- 
atures, but required an extrapolation to avoid finite size 
effects. In that work, chains of size up to iV = 14 were 
exactly diagonalized, so if we improve the linear algebra 
routines used in our implementation it should certainly 
be possible to do QBP with an Iq ^ 14. It may be pos- 
sible to do QBP with lo ~ 30 using DMRG techniques 
to avoid keeping all of the eigenstates of H and instead 
truncate to some smaller subset of eigenstates of H. Ex- 
trapolation of our results suggests that this should permit 
access to temperatures T/J^ 1/36. 

Transfer matrix DMRG[i^ is related to QBP proce- 
dure in that both procedures look for the largest eigen- 
value of a transfer matrix. In the case of transfer ma- 
trix DMRG, the transfer matrix comes from a Trotter 
approximation. In QBP, the transfer matrix is given 
by Eq. p2|l . Accurate results were obtained at much 



lower temperatures than here in [1^, but significantly 
larger matrices were diagonalized in that study, and the 
higher temperature results for peak susceptibility and 
specific heat do not appear to be as accurate. It is likely 
that the higher accuracy of our method at high temper- 
ature comes from the lack of Trotter error: the error 
becomes exponentially small in lo once lo becomes of or- 
der J(3. We now describe a procedure that combines 
some of the ideas of QBP and transfer matrix DMRG. 



Introduce M copies of the system, each with density ma- 
trix exp(— /3/i(^)/M), so that the joint density matrix 
is exp(— /3/iW/M) ® ... ® exp(-/3/iW/Af). Let Pi be 
the operator that cyclically permutes the value of the 
spin on site i between the M different copies. Then, 
tv{PiP2...PN exp(-/3/iW /M) (g) ... (8)exp(-/?/iW/M)) = 
trexp(— /S/i'^^)). Let p be a 2'"*^ dimensional matrix. 
We will define a QBP procedure such that after mt it- 
erations, p is approximately proportional to the reduced 
density matrix tri...„., {Pi...Pn„ exp(-/3/i('o-i+"« V^) 
... exp(-/3/i('«-i+"")/M)). Define O by Eq. (HID, for 
H — at inverse temperature (3/M, and map 

p Z-^T'^i^ tri(Pi(0 0) ... ® 0)p(Ot ® ... ® O^)) ® 

]i,„+i)r. (16) 

For M = 1, this reduces to the QBP implementation 
described here, while for lo = 2, this becomes very similar 
to the transfer matrix used in transfer matrix DMRG. 
The question is whether for M > 1, Zq > 2 more accurate 
results can be obtained, possibly using DMRG to find 
the fixed point p of this transfer matrix. 

QBP can be directly applied to finite size chains and 
to infinite or finite trees. The ability to compute real- 
space correlation functions and handle translationally 
non-invariant systems are advantages of this method, and 
in future this method will be applied to disordered sys- 
tems where TDMRG will have problems. Probably the 
most interesting question is the possible application of 
QBP to higher dimensional systems, by replacing the 
higher dimensional lattice with a Cayley tree or Husimi 
cactus with the correct local structure. 
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